Exact theory of dense amorphous hard spheres in high dimension 

I. The free energy 
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We consider the theory of the glass transition and jamming of hard spheres in the large space 
dimension limit. Previous investigations were based on the assumption that the probability distri- 
bution within a "cage" is Gaussian, which is not fully consistent with numerical results. Here we 
perform a replica calculation without making any assumption on the cage shape. We show that 
thermodynamic functions turn out to be exact within the Gaussian ansatz - provided one allows 
for arbitrary replica symmetry breaking — and indeed agree well with numerical results. The ac- 
tual structure function (the so-called non-ergodic parameter) is not Gaussian, an apparent paradox 
which we discuss. In this paper we focus on the free energy, future papers will present the results 
for the structure functions and a detailed comparison with numerical results. 

I. INTRODUCTION 

Hard spheres in the limit of large spatial dimensions provide us with an opportunity for an analytic solution covering 
many aspects of the liquid and glass physics [1-4]. The reason why this limit is solvable is geometric: consider three 
spheres A, B and C, with AB and BC in contact, respectively. What are the chances that A and C will themselves 
also be in contact? In high dimensions, vanishingly small. This led to the realization [1] that all terms in the virial 
expansion above the second could be neglected in high dimensions, as they involve geometrically heavily suppressed 
"coincidences" , leaving one with only the two first terms of the series for the entropy 

S[p(x)] = J d d xp{x)[l - \ogp{x)] + i J d d xd d yp(x)p(y)f(x - y) . (1) 

Here f(x) = e~ v ^ - 1 = -0(D - \ x\) is the Mayer function of the Hard Sphere potential v(x), which is infinite for 
|ar| < D and zero otherwise. D is therefore the sphere diameter. For the liquid phase, one has p(x) = p, a constant, 
and the above expression gives the liquid entropy; non-uniform phases are described by solutions of the stationarity 
equation dS/dp(x) = that are not constant 1 . The liquid phase stays mctastable at higher densities when one 
expects the thermodynamics to be dominated by a modulated, crystalline phase, which is however only known in 
small dimensions [5]. 

Even neglecting the crystalline phase, one expects that at some density, there is the possibility of a thermodynamic 
glass transition into a phase where one or several spatially-dependent non periodic solutions dominate. The conceptual 
(and practical) method to neglect the crystalline phase and uncover a possible liquid-glass transition was proposed 
years ago [6, 7]: one studies the system perturbed by a spatially random external field - whose function is to kill 
the crystal and select one of many amorphous solutions. One in fact computes the average over the "pinning" field 
realizations, and then continues the solution to zero field intensity. This program has been followed for the hard 
sphere case [8] (see [9] for the state of the art). The inclusion of a random field brings about the problem of treating 
quenched averages, and this has been done using replica methods (i.e. introducing m copies of the same system). One 
ends up with a truncated virial expansion [9] 

S[p{x)] = J dxp(x)[l - log />(£)] + i J dxdyp{x)p(y)f(x - y) (2) 

where now the density p(x) is function of m coordinates in d dimensions x = {x a } = {x±, ...,x m }, and one has then 
to analytically continue over m from integer values of m to non-integer ones. The stationary points of Eq. 2 satisfies 



1 This statement is exact at low densities. At higher densities Eq. (1) would have corrections, however these corrections are exponentially 
small in the density region we consider in this paper [3, 4]. 
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the equation 

log /3(a)] = / dyp(y)f{x-y) , (3) 



This equation has also a probabilistic meaning, as discussed in [10]. Indeed the same equation appears in the study 
of a hard sphere model on the Bethc lattice [11]. 

Even though equation (3) is (morally) exact in large dimensions, in order for it to be useful we need an expression 
for p that contains m explicitly and allows to continue the results for real m. An approach to do this is to propose a 
"Gaussian ansatz" for p(x) [8, 9]: 

pm -d/2 f 1 1,m \ 

= cxp [- w !>■ - Xb) ) (4) 

and to extremize (2) with respect to the parameters in (4). Having to content oneself with the Gaussian ansatz might 
seem somewhat disappointing: we have payed the price of going to unphysically high dimensions in order to have an 
exact answer, and now we do not even have this. 

The purpose of this paper is that of obtaining the exact solution of (3) in the high-dimensional limit without 
assumptions. We find that the Gaussian ansatz turns out to be, in a sense, exact: it gives for large dimensions the 
exact values for the thermodynamic quantities. The reason, whose consequences we shall develop below, can be seen 
as follows: a generic replica problem is written in terms of the order parameter p(x), or, alternatively, of the tensors 
(x a Xb),{x a XbX c ), (x a XbX c Xd), ■•■ The solution involves an ansatz for replica tensors of all degrees, which makes the 
problem analytically hard. However, in our case, the x a are vectors in d-dimensional space, and we are looking for 
a solution that is statistically translationally invariant and isotropic. The only possibility with these properties is 
that p(x) depends exclusively on the \x a — Xb\ 2 = q aa + Qbb — 2q a b, where we have introduced the scalar products 
Qab = x a ■ Xb- All e?-dimcnsional integrals may be expressed as low dimensional integrals in terms of the q a b, with 
volume factors, a simple generalization of spherical coordinates. As we shall sec, the dimensionality appears explicitly, 
and the limit of large d may be taken in a straightforward way, by saddle point evaluation of the integrals. 

It will turn out, however, that the "cage shape" is not Gaussian, as already observed in simulations [12], but 
may be calculated exactly for large d within this framework. The fact that the Gaussian approximation gives the 
correct result for thermodynamic functions but by itself does not give the right cage shape may be understood with 
a simple example. Consider a single particle in a d-dimensional spherical potential V, which for convenience we write 
as PV = d U(\x\ 2 /d). Let us denote q = r 2 = \x\ 2 . Clearly, the exact Gibbs distribution is 

p -dU(q/d) 1 f 

p(q) = -— with AA G = i dq q^' 2 e~ dU ^ (5) 

■A/o ^ J 

The entropy S = — J dx phip is easily evaluated and gives: 

f -dU(q/d) 

S=- dq q^' 2 {-d U(q/d) - lnj\f } (6) 

For large dimension d, the integrals for Af and for S are dominated by the saddle point q* which maximizes 
hkiq — U(q/d). The entropy is then, to leading order in d: 

S = -dU(q*/d) - lnAA - -dU{q*/d) + dU(q*/d) + ^\nq* = ^\nq* (7) 

Suppose instead that we had done this calculation approximately, proposing a variational Gaussian distribution 
pa = e - ' 35 ! /( 2A ) /Mg- The same calculation as before gives for the entropy: 

f d 

S G = - dx pg^Pg ~ 2 ln?A ( 8 ) 
where qA maximizes ^mg — q/(2dA). We now have to fix the parameter A by minimizing the free energy: 

/d 
dx [-0V - In p G ] PG - -dU{q A /d) + -lnq A (9) 

The minimum is clearly attained by q A — q* , by definition of q* . All in all, we have to choose the value of A such 
that qA = q* , the "true" saddle point. Entropy and free energy give, for this value, the correct results Sq = S and 
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Fq = F. We note that the only purpose of the Gaussian ansatz at this stage is to fix the correct value of q* that 
dominates all integrals, a purely geometric feature: indeed, any other ansatz (e.g. a delta function) would have given 
the correct result by a similar argument. This is akin to the equivalence of different thermodynamic ensembles in 
infinite dimensional space. 

If in the same example we are interested in the radial "cage" function p(r), the result will be Gaussian if we use 
pa, while it will generically have tails ~ e~ du that are not Gaussian in the exact case. In other words, the tails of 
the "cage" distribution are large deviation functions with (large) parameter d. They could also be computed in the 
Gaussian approximation by evaluating the free energy as a function of the intensity h of an additional potential hx: 
as usual large deviations control the thermodynamics in presence of an external field. 

This paper is organized as follows: in section II we write the replicated Van der Waals entropy in coordinates 
corresponding to the scalar products, as described above. Sections III and IV are devoted to the calculation of the 
Jacobians of these changes of coordinates, that will play the role of the term r d ~ x of polar coordinates in the example 
above. (We need two Jacobians, corresponding to integrations in spaces of m vectors {x a }, and to 2m vectors {x a , y a } 
- as required by the Mayer function - in terms of the corresponding scalar products). In section V we compute the 
Mayer function in terms of these coordinates. Thus, we obtain a complete expression (section VI) for the entropy, in 
terms of low-dimensional integrals, with the dimension d appearing as a parameter. In sections VII and VIII we do 
the analogue of the previous paragraph: we compute the thermodynamic functions using the Gaussian ansatz and the 
exact solution taking saddle points that become exact as d — > oo. Both results coincide, thus validating the Gaussian 
ansatz. 



II. THE REPLICATED VAN DER WAALS ENTROPY 



The starting point of our calculation is the free energy of a replicated liquid, where each atom is replaced by a 
"molecule" made by one atom per each of the m replicas [8, 9]. We denote by x = {x\, ■ ■ ■ ,x m } the coordinate of 
such a molecule, each x a being a vector in <i-dimensional space. We assume that in the glass phase the molecule is well 
defined, the typical distance between atoms in a molecule being of order \fA which is small at the glass transition [9] . 
Note that this is a non-trivial assumption: molecules might dissociate, especially close to the glass transition, and 
lose their identity. However, we will show self-consistently in the end that this is not the case, at least for d — > oo. 
Taking into account this effect in finite dimensions might be a non-trivial task. 

The liquid state is described by a single copy of the system, m = 1, with uniform density p. When d — > oo, its 
entropy (per particle) is given by Eq. (1) for p(x) = p, which corresponds to the Van der Waals mean field equation: 

-a, = ^ = 1 - log, + ^ = 1 - logp + 2-" V • (10) 

Here Vd is the volume of a sphere of unit radius in d dimensions, and D is the sphere diameter, and we introduced the 
packing fraction ip = pVd(D/2) d . It has been shown in [1, 3, 4] that Eq. (10) is exact in the limit d — > oo, provided 
2 d ip does not grow exponentially with d. If this is the case, the other virial corrections are exponentially suppressed 
in d. 

In the replicated liquid, atoms within a molecule can overlap, while atoms of different replicas belonging to different 
molecules have the normal hard sphere interaction. If y/A -C D, the molecule-molecule interaction is similar to the 
normal hard sphere interaction and one can repeat the analysis of [3]. The replicated liquid with integer m > 1 can 
thus be described in terms of a replicated Van der Waals entropy given by Eq. (2) (see [9, 13]), where p(x) is the 
single molecule density, normalized to J dxp(x) = p where V is the system volume, and f{x — y) is the replicated 
Mayer function that describes the molecule-molecule interaction: 



f(x -y) = -l + Y[ 6(\x a - y a \ - D) . (11) 



£1=1 



We wish to make use of the homogeneity of the molecular liquid, which implies that we have to consider a generic 
translationally and rotationally invariant form of p{x). Let us start with translational invariance. We can perform a 
change of variables, X = m _1 x a , u a = x a — X , where X is the center of mass of the molecule and u a are the 
relative displacements with respect to X . Then 

dx = dX du m d 8 u "j = dX Vu > ( 12 ) 
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and translational invariance implies that p(x) does not depend on X . We obtain 

)J{u-v) , 

(13) 



Vup{u)[l-\ogp{u)] + - I 



f{u-v) = / dXf{X + u-v) 



where X + u means adding X to each component of u. 

Next we consider rotational invariance, which implies that p(u) and f(u) are functions of q ab — u a ■ Ub only. Let us 
define the following quantities: 

q a b = u a - u b 

Pab =v a -v b (14) 

Tab =U a -V b 

Note that translational invariance implies that JZaLi Qa,b = ^2T=i l^b = for all row and columns. Denoting dq = 
Ha<q d Q°-b, we have 

S[p{q)]/V = [Vup(q){l- log p{q)} + I [ VuVv p(q) p(p)J (q + p - f -f T ) 

f 1 f - ( 15) 

= / dqj(q)p(q)[l- log p(q)] + - / dqdpdrK(q,p,r)p(q)p(p)f(q+p-r~r T ) . 

Here we introduced two Jacobians J{q) and K(q,p,f) that describe the change of variables from u to q, and from 
(u,v) to (q,p,f) respectively. We compute them in the next sections. Before doing that, note that translational and 
rotational invariances imply that the integrals are reduced from ~ md variables u to ~ m(m— l)/2 variables q. This 
is crucial because now the number of integration variables docs not grow with d and we can use saddle-point methods 
when d — > oo. 

III. THE JACOBIAN J 
A. Definition 

The Jacobian J(q) is defined as: 

l,m _ / ra \ l,m 



Vu Y[ S(q ab - u a ■ u b ) =m d du S I u a J ]~[ <5(q a6 - w a • u 6 ) 

a<6 ^ \a = l / o<b 

m / m \ r l,m— 1 

= TO d (5 (? ab dm--- dum-i \\ 5(q ab - w Q • u 6 ) , 

1 \ L 1 / J _^-L 



(16) 



where the second line is obtained easily by manipulating the delta functions. 

The delta functions take into account translational invariance. The last term instead takes into account rotational 
invariance, and it can be shown that 



l,m— 1 



( dm--- du m . t f[ S(q ab - u a ■ u b ) = C m , d e ^-™)i°gdct^- ; (17) 

a<b 

where q a > b is the (m — 1) x (m — 1) matrix that is obtained by removing from q the a-th row and the &-th column. In 
fact, Eq. (17) can be thought as a change of variables from a d x (m — 1) matrix U = {u\, • • • , u m -±} whose columns 
are the coordinates of the vectors u\ ■ ■ ■ u m -i, to the (m — 1) x (m — 1) matrix q m < m — JJ T (J. The corresponding 
Jacobian is given by Eq. (17), see [14]. 
Using this, we get the final result: 



in / m 



J(q) = m d C mtd Y[6[Y,<lab) e ^ d - m ^ dct ^ m . (18) 



a=l \b=l 
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Note that this form of J is consistent with the fact that the choice of the m-th row and column in (16) is arbitrary: we 
could have chosen any other row and column. But because the matrix q has the property that the rows and columns 
add up to zero, it has the property that dct q a ' b = det q m ' m for any a, b. 
The normalization constant is 

d 

C m ,d = 2 1 ~ m Ofc , 

k=d-m+2 Qg'j 

2?r d / 2 

^ d = t7~j /o\ ' ^ ae ^-dimensional solid angle) 

1 (a/ 2) 

as we show in next section. 



B. Calculation of the normalization constant 



Here we compute the normalization constant C mtd . First we note that one can compute it directly for m = 2 and 
m = 3, using polar and bi-polar coordinates respectively: 

r _ V d 

L'2,d — — , 

il d \l d -l 

c 3 , d - — - A — • 

This already hints strongly at the form (19). 

Next we perform an asymptotic computation for large d at fixed m. For this, we write (dropping for convenience 
the superscript (m, m) on the matrix q that here we consider here to be a generic (m — 1) x (m — 1) matrix): 



l,m-l 



( 27r )d(™-i)/2 = I dui ... du m - ie -i i 1 u "- u " 

= / dq / dui---du m -i e - '^^ 1 qaa S{q a b-u a -u b ) (21) 

J ^ a<b 
= C m>d J dq e i («i-*n) log de* «- f S^l 1 qaa 

The latter integral can be evaluated by a saddle point. Using -^-^ log dct q = (q^ 1 )^- the stationary equation reads: 

-1+ (d- mjq- 1 = , q = (d-m)I. (22) 

Next we expand 

q=(d-m)I + i (23) 
and we obtain, truncating the expansion in t at quadratic order: 

(27T) d(m " 1)/2 = C md J dt e 5( d -™){(" i - 1 ) lo S( d -"0+Tr[t/(d-m)-tV2/(d-m) 2 + ---]}-iEr=Y t aa-5( d -"0(™-l) 

= C m4 {d - m )(d-m)( m -l)/2 e -i(d-m)( m -l) J di g-jp^y & (24) 

- C m , d (d - m )(^)(rn-l)/2 e -Ud-rn)( m -l) 2 ( m -l)/ 2^^ _ m) " l(m - 1)/2 

We therefore get 

C m d = 2-i( m - 1 '( 2+m " 2d )e5( d - m )( m - 1 ) 7 r"3 (m ~ 1)(Tn ~ 2d) (d - m,)!™^- 1 )"! ^C" 1 - 1 ) (25) 
It is easy to show that the limit for d — > oo of this expression divided by Eq. (19) is given by 1. 
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IV. THE JACOBIAN K 
A. Definition 



The Jacobian K(q,p,f) is defined as: 

/ l,7n l,m l.rn 

VuDv Y[ 5(q ab - u a ■ Ub) 5(p ab - v a ■ v b ) Y[ 5(r ab ~ u a ■ v b ) 

a<b a<b a,b 

r ( m \ ( m \ 1,m 



m 2d 



(26) 



Va— 1 / \a— 1 / a<b a<b a,o 

m / m \ m / m \ / 1 , m — 1 \ m — 1 / m \ m— 1 / m 

= m 2d n<5 £ to n<5 $>* W r TOm - £ ^ n <m i> n * x><* ) 

a=l \fc=l / a=l \6=1 / V a,b J a=l \6=1 / 6=1 \a=l / 

/l,m— 1 l,m — 1 l,m — 1 

dui • • • du m _idvi ■ ■ ■ dv m -i _ <5(g a6 - u a ■ u b ) S(p ab - v a ■ v b ) S(r ab - u a ■ v b ) 

a<b a<b a,b 

where the last line is obtained easily by manipulating the delta functions. 

We can define again a matrix U = {iti, • • • , u m _i, Vi, ■ ■ ■ , v m -i} of size d x 2(m — 1) and a matrix Q = U T U of 
size 2(m — 1) X 2(m — 1), such that 



Q = 



%rn ,m\T , m 



(27) 



is obtained from the matrices q,p,r from which the m-th line and column has been removed. 
Clearly we can write, using Eq. (17) and calling U a the columns of the matrix U: 



l,m — 1 l,m — 1 



/_L ■ r r ' _L _l_ ■ I I L _l_ J. » ( lb A 

dm ■ ■ ■ du m -idvi ■ ■ ■ dv m -i _ S(q ab - u a ■ u b ) _ _ 5{p ab - v a ■ v b ) _ _ 5{r 
ab u a ' ^b) 

a<b a<b a,b 

l,2(m-l) 

= / dUi ■ ■ ■ dU 2 ( m -l) II S ^ ab ~ Ua ' Ub ) = C 2m-l,d 



(28) 



;(d-2m+l) log dot Q 

a<b 

Using this, we get the final result: 

K(q,p,f) =m 2d C 2m . hd e i(«*-am+l)logdet<5 

m / m \ rn / m \ / 1 , m— 1 \ m — 1 / m \ m — 1 / m ^ 

x n s ( E ?a6 ) n ^ ( E po6 ) ^ ~ e n ^ e rafc ) n * e ^ 

a=l \fc=l / a=l \b=l / \ a,b ) a=l \6=1 / 6=1 \a=l / 



V. THE REPLICATED MAYER FUNCTION / 
A. General expression 

We now investigate the replicated Mayer function f(u), which is defined as 

7(u) = J dxS^-l + f[eQX + u a \-D)^ = - J dXe(D-mm\X + u a \) 



(29) 



(30) 



The u a are m vectors in d dimensions. In the following we assume that d > m. A remark that will be useful in the 
following is that when all u a = 0, / = —Va\D d , while when the distance between each pair \u a — uj,\ > D, we have 
/ = -mV d D d . 
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We define X\\ as the part of X that lies in the hyperplanc defined by the u a and X± the orthogonal part. Then, 
recalling that is the d-dimcnsional solid angle and Vd = £ld/d, 

J(u) = - dX 6{D- min \X + u a \) = - dX 9(D 2 - min \X + u a \ 2 ) 
J a J a 

= -Jd m X l{ d d - m X ± 0(D 2 -nun{|X|| + u a \ 2 + \X ± \ 2 }) 

/poo 
d m X\\ J dxx 4 -™- 1 9{D 2 - x 2 - min \X\\ + u a \ 2 ) 



-0 d _ m / d m Xn I dxX d - m - X 



(31) 







= -Vd- m [ d m X\\ (D 2 - min + u a \ 2 )( d - m V 2 0(D 2 - min |X„ + u a \ 2 ) 

J a a 

= -V d - m j drX\\ e d - m (D 2 - min|X|| +u a \ 2 ) 



where we defined the function 



Qd-m(x) = x^ d - m ^ 2 0(x) . (32) 



While the above formula is always valid as long as d > m, for large d the last integral is dominated by the points 
where min a \Xu + u a \ = 0, which means that = — u a for some a (each value of a defines a different saddle point). 
Observing that / = — VdD d when u a = Va, we can also write: 



. / <f"Xj| e d ~ m (D 2 - min a + u a \ 2 ) 



B. Evaluation of / for d — > oc 

Let's consider first the case where the vectors u a are very large. In this case, if we write Xii = — u a + e, for small e 
the minimum min a \Xn + u a \ will still be assumed in the same value of a than for e = 0, hence min a \X\\ + u a \ = \e\. 
Then we get 

m « 

l{u) ~ Vd— jtj 5^ / <^ - \e\ 2 ) (d - m)/2 6(D - \e\) 

a=l ^ 

~ -mVd- m ^ra f dee m -\D 2 - e 2 r-^' 2 ( 34 ) 

dT(d/2) 

which implies that f(u) is a constant exactly equal to minus the volume of m hyperspheres, — m x VdD d . Note that 
the integral over e is dominated by a saddle-point at e ~ l/v 7 ^, as it can be easily checked. On the other hand, in 
the limit |it a | = for all a, we trivially obtain f{u) = —V d D d , the volume of one hypersphere. 

Therefore, the region where / has a non-trivial dependence on the u a is where the u a have a length proportional 
to l/Vd. We can define u a — x a D/\Jd — m and Xh = eD/y/d — m, and we can write from Eq. (33): 



f(u) = ~V d D d 



f ( jm e ( I _ min a |e+:c a | 2 \ ^ m '' 2 Q ( <\ _ min a |e+:c a | s 
J I d—m ) \ d—m 



-ViD 



d J d m e milla l e + x »l 2 



(d-m)/2 

(35) 



11^12 



M 2 



which is still of the order of VdD d times a non-exponential factor that depends on the u a . 



We therefore conclude that f(u) has the following scaling form when d^oo: 



7(u) = -V d D d T hHI^L u 



f Jm, „— A min a U+rrJ 2 r Jm 



where 

Note that when all x a — 0, J- — 1 as it should; and when each distance \x a — x b \ — > oo, J 7 — > m. 

VI. ROTATION ALLY AND TRANSLATIONALLY INVARIANT EXPRESSION OF THE REPLICATED 

VAN DER WAALS ENTROPY 

A. Exact expression of the entropy 

Before proceeding, let us collect here the results obtained up to this point. We wrote the replicated Van der Waals 
entropy, taking into account explicitly rotational and translational invariancc, as follows: 

SW)\/V= f dqJ(q)p(q)[l-logp(q)} + ^ f dqdpdrK(q,p,r)p(q)p(p)J(q + p - f - f T ) (38) 



where 

m / rn 

J{q)=m d C m4 X[5 $>ab) e ^ d - m ^ dot « m - m (39) 

a=l \6=1 / 

(here q m ' m is obtained by removing the m-th line and column from q) and 

K(q,p,f) =m 2d C 2m _ M e l(d-2 m+ i)iogdet<3 



m / rn \ rn. / rn \ / l,m— 1 \ m — 1 / m \ m— 1 / m ^ 

x IT 5 qab ) IT 3 Pab ) a Tmm ~ J2 rab ) Yi 5 \ J2 rab ) Yi 5 [J2 rab 

a=l \6=1 / o=l \6=1 / \ a,b I a=l \b=l / 6=1 \a=l / 



(40) 



(here Q is obtained from the matrix 



q r 

~T - 

r p 



by removing the m-th and 2m-th lines and columns) and 



C md — 2 1 ~ m Y[ ^fc ~ 2~^ m ^ 1 ^ 2+m ~ 2d ^e^ d ~ m ^ m ~ 1 ^~^ m ~ 1 ^ m ~ 2d \d — m)^"^" 1-1 ) - ^" 1 - 1 ) (41) 



k— d— m+2 

B. Equation for p(q) 

Remember that p{u) is normalized by V~ x J dxp(x) = J T>up(u) = p. Starting from Eq. (38) and differentiating 
with respect to p(q), adding a Lagrange multiplier to ensure normalization, we obtain the equation: 

J{q) logp(<?) = J drd P K(q,p, r)p{p)J(q+p- f - r T ) + XJ(q) (42) 
Recall now that by definition: 

l.m l.m 



/ dfdpK(q,p, f)p(p) = I dpdf I VuVv TJ S(q ab - u a ■ u b ) TJ S(p ab - v a ■ v b ) TJ 5{r ab - u a ■ v b ) p(y) 

1,771 

= / VuDv Y[ 8{q ab - u a ■ u b ) p(v) = pj(q) 
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Therefore we can write the equation for p(q) as follows: 

, / ~\ , , JdrdpK(q,p,f)p(p)J(q + p-r-f T ) 

J drdpK{q,p,r)p(p) 

where obviously the delta functions involving q in the expression of K have to be formally simplified between numerator 
and denominator. The multiplier A is determined by the normalization condition. 



VII. THE GAUSSIAN ANSATZ 



Before moving to the general case, we show here that the computation above gives back exactly the results of [9] if 
a Gaussian ansatz is made for p(q). Using this Gaussian ansatz, we evaluate Eq. (38) using the saddle point method. 
In the next section we will show how this result can be obtained in fully generality. 

We observe that because neither log p(q) nor f(q) arc exponential in d, the saddle point is only determined by the 
Jacobians and by p(q) in both terms of Eq. (38). Therefore, f = at the saddle point, and q — p = q sp , where q sp is 
the point where the exponential factor in J(q)p(q) is maximum. Substituting this saddle point in Eq. (38), we obtain 

S[p(q)]/N ~ 1 -logpOT) + ^f(2q s P) 

~ 1 - log p(q s P) - 2 d -\T (^2^ 

where the second line is obtained from the first by using Eq. (36). 



A. The Gaussian form of p(u) and the entropic term 



The Gaussian ansatz has the following form: 

P( ' (2nA)( m - 1 W 2 P[q> (2 7 r J 4)( m - 1 ) rf / 2 1 ' 

The first task is to compute the saddle point value of q that dominates all the integrals. We have, using the delta 
functions contained in J(q) to manipulate the exponential term in p(q): 

/„ m /to \ 

dqJ{q)p{q) (X / dq JJ S ^ q ab \ e h(d-m) log detq^-^^S, 1 laa+llT' 1 ( 47 ) 

At this point the integral over q am is eliminated by the delta functions and we are left with the (m — 1) x (m — 1) 
matrix q m ' m . The saddle point equation is, for d — > 00 and a, & =!,-■■, m — 1: 



The matrix is easily inverted and we obtain 



(r 1 )Z = ^ I (l + 6 ab ). (48) 



q^ b =dA{5 ab --) (49) 



It is easy to show using the conditions X^fcli Qab = imposed by the delta function that the formula above holds for 
a,b — 1, • • • ,m. Indeed the saddle point values satisfy = (u a ■ Ub) (where the average is over p(u)) so the same 
result could be obtained from a direct computation. We get 

, c „n , d, (m — l)d (m—l)d, , , 

1-logpOT = l-logp+-logm+ V 2 ' + l - ' log(2nA) (50) 

which is the same result that can be obtained by an exact computation, the integrals being Gaussian in this case [9] . 
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B. The interaction term 

Next, we compute the term involving / in the saddle point. Let us start with the following observation. Because 
f(u) depends only on q thanks to rotational invariance, all values of u that correspond to the same q give the same 
value of f(u). This means that if we want to compute /(<?), we can do that by choosing our favorite configuration of 
u that corresponds to the chosen q. 

Therefore, for the saddle point (49), we can choose any u sp that satisfy the conditions 



E< p = 



a=l 



sp 
lab 



Snh - ~ 



(51) 



Remember that u sp are d-dimensional vectors. A good choice is {u^ p ) b = \J Ad (5 a b — tJt) for their first m components, 
b = 1, • • • , to, and zero for all the other components. 

We therefore use this configurations of the u to compute f(u). We further define A by A = D 2 A/d 2 . Therefore, 
the corresponding variables x that appear as the arguments of Eq. (37) in the saddle-point Eq. (45) have to satisfy 



rr S P . rJ. S P 

x a x 6 



-,2q2 ~2A\5 ab -- 



D 



in 



and therefore 

With this choice, a short computation shows that 



2A S, 



niin| e + *-| 2 = 5> b ) 2 -^£ 

n. — * rrj ^ — * 



6=1 



6=1 



e + 2V2Amine a + 2A ( 1 

TO 



Therefore 



J"(5 sp ) 



to e 



2tt 



e 2 



min Q |e+< p | 2 



dc 



27T 



e 2 



drj 



e 2 ' 



p oo 


de 1,2 




—=e 2 


J — oo 


V2^ 








+ erf 







'2A< 



2tt 



1 1 + erf V2A - £m 

2 I I v^to 



_j V2tt 
It is useful to define 

e? m (l) = i-jr(x sp ) = i-m 



'2A- 



de _i c 2 
=e 2 



2tt 



5 1 + -0^T 



m— 1 



(52) 



(53) 



(54) 



(55) 



(56) 



C. The Gaussian result 



The final result of the Gaussian computation is therefore, at the leading order for d — >• oo: 

„ r , , d, (m—l)d (m — l)d , \2-kD' 2 A\ , « r „ , ~. n . . 

S[p(g)]/JY = l-bgp+-logm + V 2 + 2 log I — j -2 d -V[l-gm(A)] . (57) 

It coincides exactly with the result of [9]; the explicit expression of Q m (A) that was given in [9] is different, but it is 
exactly equivalent to the present one (actually the present one is much easier to compute numerically). 
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VIII. THE GENERIC 1-STEP REPLICA SYMMETRY BROKEN SADDLE POINT 

Finally we analyze the structure of the generic solution for p(q) in the large d limit without assuming a Gaussian 
form for the density function. Like in the Gaussian case, we want to evaluate the integrals in Eq. (38) and Eq. (44) 
via a saddle point. We will show that we will recover the results coming from the Gaussian case. 

A. Structure of the saddle point 

First we have to derive the saddle point equations. Let us suppose that 
p(q) = e" n «) 

/ m— 1 m— 1 1,TO— 1 

U)(q m ' m ) = il q ab , q am = - ^ q ab , q mb = - ^ q ab , q mm = ^ q a , 



(58) 



b— 1 a— 1 a.b 



The 1-step replica symmetric breaking (1RSB) solution consists in assuming that Q(q) has a replica symmetric (RS) 
structure. Indeed, this corresponds to 1RSB because the present real replica scheme describes what happens inside 
one of the 1RSB blocks [7-9]. 

First we want to determine the saddle-point value of q, that dominates the normalization of p(q). We have 

p= I dqJ(q)p(q) oc / dqJ[S I X> 6 e^-^ 1 ^^^-^"^ . (59) 

a=l \6=1 / 

Therefore the delta functions allow to eliminate the m-th line of q ab , and the saddle point equations for the remaining 
variables q ab , with a, b = 1 • • • m — 1, are determined by the maximization of the exponential factor for d — > oo: 



,~-i,sp du{q) dfl dfl dfl dtt 

) ab = —j = 3 + -j -j -j ■ ( 60 ) 

1 aq a b aq a b aq mm aq am aq mb 



Because 0(g) has a RS structure, we have = Slo(<z) + 2n 1 (q) ^ abj anc ^ 

d \ m 



(61) 



Consider now the integral: 



p = I dqdpdfK(q,p,f)p(q)p{p) . (62) 



A very similar procedure leads to the following saddle point equation: 



d - i d r 



2^ 2 



o 2^(1 + <U) 



(63) 



whose solution is r ab = and q ab = p ab = q'^, both equal to the solution of Eq. (61). Finally, we are interested in the 
integral entering in Eq. (44), that is the same as the last one but without the condition on q. One obtains the same 
as Eq. (63) but without the upper left block of the matrix. However, the equation for f is still solved by f = 0, then 
the equation on p decouples from q and leads to the same solution as in Eq. (61). 

We conclude that the integrals in Eq. (38) and Eq. (44) can be evaluated via a saddle point and the result is 

Slp(q)]/N~l~\ogp(q s n + pm SP ) 

( d \ ( 64 ) 

~\-\ogp{q*v)-2 d -\T{—{lq s v 
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and 

logp(q) = X + pf(q + q sp ) 



A-2<V^(^(<7 + <f P )) 



d . „,„A (65) 



Therefore, to conclude the calculation we have to determine q sp and A. 

Before proceeding, two remarks are in order. First of all, when the distance between atoms in a molecule is large, 
the q a b are large, and J->mas discussed in Section V. Hence in this limit p(q) ~ exp(— 2 d ipm), and because 2 d ip oc d 
at the glass transition [9], we see that p(q) goes to a very small constant that vanishes exponentially with d. Hence, 
in the d — > oo limit the molecules are well defined in the glass phase, while for finite d there is an exponentially small 
probability of dissociation. This guarantees that the molecular liquid is a good description of the glass phase for 
large d. The second remark is that the choice of a given RSB ansatz is self-consistent. We assumed at the beginning 
a RS structure for fl(q); then we obtained that q sp has a RS structure as given in Eq. (61); and finally that, self- 
consistently, Q(q) = — log p(q) has a RS structure as given by Eq. (65). We could consider a 1RSB structure for Q(q) 
(corresponding to a 2RSB computation in the real replica scheme) and we would have obtained self-consistently the 
same structure for q sp . Saddle points characterized by many steps of RSB could be needed to describe the metastable 
states of lower density [15]. 



B. Saddle point equation 



We now have to solve the saddle point equation Eq. (61). We note that, defining A sp = fli(q sp )/ D 2 , we can rewrite 
Eq. (61) as a closed equation for the scalar parameter A sp : 



sp 
<U = 



D 2 A sp 



Sab ) , 

m . 



Furthermore, if we define a function 



then wc have 



A sp = — n x 

D 2 



h(A) = n 



D 2 A sp 



Sab 



D 2 A 



Sab 

m 



^M sp ) = E ^ sp ) ? ('* --) =E( 

dA ^ dq ab d \ mj 



ab 



20j (qw) 



D' 1 



Sab — T \5ab = 



1 \ _ dD 2 (m - 1) 



2fii (q sp ) 



therefore 



and the equation for A sp becomes 



A sp 



dD 2 (m- 1) 
2h'(A sp ) 

d(m - 1) 



2h'{A sp ) 

The function h(A) can be computed by using Eq. (65), that gives Sl(q) = -\ + 2 d ipT{^(q + q sp )). Then 

h(A) = -X + 2 d v F 



(66) 



(A + A sp ) [S ab -- 
m 



(67) 



(68) 



(69) 



(70) 



(71) 



The computation of this function can be done exactly as we did in the Gaussian case, in section VII B, and leads to 
the following result: 



(A + A SP ) ( Sab - - 

m 



= l-9 m [{A + A sp )/2] . 



(72) 
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Therefore 

h'(A sp ) = -2 d -\g' m {A s P) , (73) 

and finally the equation for A sp is 

A sp = i ^— . (74) 

2*<pQ' m (Asp) 

It is easy to check that this equation is exactly the same that is obtained by maximization of the Gaussian free entropy 
Eq. (57). Therefore, the generic saddle point equation coincides with the Gaussian one. 



C. The computation of A 



The last ingredient that we need to compute the replicated free entropy in the generic case is the value of A. Indeed, 
combining Eq. (64) and Eq. (65), and recalling that JF (j^2q sp ) = 1 - G m {A sp ), we obtain 

S[p(q)]/N = 1 - A + 2 d "V [1 - G m (A sp )} . (75) 
The factor A has to be computed by imposing the normalization of p(q): 



/„ m / m 

dqJ(q)p(q) = e x m d C m , d / dq JJ I H ]T] qab 

a=l \b=l 



i(d-m) log dctq m < m -2<V^(^-(<?+<f 3p )) 



(76) 



It is enough to evaluate the integral at the saddle point level to get the part of A that is proportional to d. Corrections 
to the saddle-point give corrections to A that are at most proportional to \ogd and will be neglected here. We get 
from Eq. (66) 

logdct(g sp ) m ' m = (m - 1) log(D 2 A sp /d) - logm . (77) 
Recalling that from Eq. (25) we have, neglecting corrections proportional to log d: 

logC m , d = ^(m-l)log(2vre) - |(m- l)logd , (78) 

and once again that T [^2q sp ) = 1 — G m (A sp ), we obtain the equation for A: 

log p = A + dlogm + |(m - 1) log(2^e) - ^(m - 1) log d + | [(m - 1) log(D 2 A sp /d) - logm] - 2<V[l - G m (A sp )] 
= A + d - logm + l(m - 1) + |(m - 1) log(2nD 2 A sp /d 2 ) - 2 d <p[l - G m (A sp )] , 

(79) 

where log d terms are neglected. Plugging this in Eq. (75) we finally obtain 

S[p(q)}/N = 1 - logp+ ~ logm+ ~(m - 1) + |(m - 1) \og(2irD 2 A sp / d 2 ) - 2 d ~\ [1 - G m (A sp )\ , (80) 
which coincides exactly with the Gaussian result, Eq. (57). 



IX. CONCLUSIONS 



We have derived the large dimensional limit for the statistical mechanics of dense amorphous hard spheres. We 
have shown that a so-called Gaussian ansatz gives the correct result for the thermodynamic functions. The results 
previously obtained with such an ansatz are thus validated, at least near the transition where a one-step replica 
symmetry breaking scheme is expected to suffice. We notice that the same computation would also apply to the 
Bcthc lattice model of [11] in the high coordination and large dimension limit. 
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There is however reason to suspect that the actual exact result has an infinite number of breakings, at least in the 
limit of high pressure: first of all, the generic situation with systems which have a transition to a one-step solution 
is to have a further transition to a phase with more - eventually infinite levels of replica symmetry breaking. More 
physically, we know that hard spheres at large pressure develop many soft vibrational modes [16-18] due to isostaticity. 
This is true not only of the equilibrium states, but also of the metastable "J-point" states. Now, 1RSB equilibrium 
states have a spectrum with no soft modes - and this is true of all but the very highest metastable states. A full replica 
symmetry breaking scheme would naturally bring in soft modes, as happens for example in the case of spin-glasses. 
Perhaps the transition into such a phase also brings in isostaticity at high pressure, something the 1RSB solution 
displays only for the equilibrium states [9]. The study of the stability of the 1RSB solution is under investigation. 

As mentioned in the simple example of the introduction, the large d calculation we have presented cannot be 
expected to yield the exact result for the cage distribution, especially its tails. A more detailed calculation, always 
within this framework and based on Eq. (65), is possible for the tails of exponentially small probability, both at the 
glass transition and at jamming. Hopefully such calculation will be able to reproduce the numerical results of [12], 
where a large non-Gaussian tail has been detected in the self part of the van Hove function, which coincides with the 
cage distribution. It was found that this tail is not reduced on increasing dimension and seems to persist even for 
d — > oo, suggesting that it could be described by mean field theory This will be the subject of a future paper. 

To conclude, let us mention that it would be nice to reproduce the results obtained here without using replicas, 
i.e. by finding directly the amorphous solutions of Eq. (1). This approach (which is also called Density Functional 
Theory or DFT) has been pursued in [2], under the assumptions that (i) the density field p(x) is the sum of Gaussians 
centered around amorphous reference positions Ri and (ii) the structure factor of the Ri is the same as those of the 
liquid. The results of [2] are close but not exactly equivalent to the ones of replica theory. It is likely that hypothesis 
(i) is not needed thanks to the same mechanism that was exposed in the introductions and is at work in the replica 
calculation. However, it is less easy to refrain from using hypothesis (ii) in DFT, and it is likely that this hypothesis 
is false. In fact, the liquid structure becomes akin to the one of the ideal gas when d — > oo, while it is likely that 
the Ri remain correlated, at least for neighboring particles (as it is clear for instance in the jamming limit of infinite 
pressure). The replica method is indeed designed to integrate over the unknown Ri to get rid of them and avoid 
hypothesis (ii). Reconciling DFT with the replica method requires a way to solve for the Ri, which has yet to be 
found. 
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